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Abstract 

The least-squares finite element method(LSFEM) based on the velocity-pressure- 
vorticity formulation is applied to large-scale/three-dimensional steady incompressible 
Navier-Stokes problems. This method can accommodate equal-order interpolations, and 
results in sy mm etric, positive definite algebraic system which can be solved effectively by 
simple iterative methods. The first-order velocity-Bernoulh function-vorticity formulation 
for incompressible viscous flows is also tested. For three-dimensional cases, an additional 
compatibility equation, i.e., the divergence of vorticity vector should be zero, is included 
to make the first-order system elliptic. The simple substitution or the Newton’s method 
is employed to linearize the partial differential equations, the LSFEM is used to obtain 
discretized equations, and the system of algebraic equations is solved using the Jacobi pre- 
conditioned conjugate gradient method which avoids formation of either element or global 
matrices (matrix-free) to achieve high efficiency. To show the validity of this method for 
large-scale computation, we give numerical results for 2D driven cavity problem at Re = 
10000 with 408 x 400 bilinear elements. The flow in a 3D cavity is calculated at Re = 100, 
400, and 1,000 with 50 x 50 x 50 trilinear elements. The Taylor-Gortler-like vortices are 
observed for Re = 1,000. 
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1. Introduction 


Although significant progress has been made in the finite element solution of incom- 
pressible viscous flow problems, development of more efficient methods is still needed 
before large-scale computation of 3D problems becomes feasible. This paper presents such 
a development. 

The most popular finite element method for the solution of incompressible Navier- 
Stokes equations is the classic Galerkin mixed method based on the velocity-pressure for- 
mulation, e.g., see Refs. [1-10]. One of notorious difficulties in the mixed method is the 
satisfaction of the Ladyzhenskaya-Babuska-Brezzi (LBB) condition which requires the use 
of different elements to interpolate the velocity and the pressure in order to obtain a sta- 
ble scheme. For two-dimensional problems quite a few convergent pairs of velocity and 
pressure elements have been developed, however most of these combinations employ ba- 
sis functions that are not convenient to implement. For three-dimensional problems, this 
difficulty becomes more severe and only rather elaborate constructions can pass the LBB 
test. Another difficulty is due to the lack of symmetry and positive definiteness of the 
linear equations arising from the classic mixed method. Iterative methods for solving this 
kind of linear system have been hard to come by. As Gresho, Lee and Sani[ll] pointed out 
“there is currently no iterative solution method which can be guaranteed to work on the 
discretized, primitive-variable incompressible NS equations”. Therefore, direct Gaussian 
elimination has been considered the only viable method for solving these systems. Un- 
fortunately, for three-dimensional problems the computer resources required by a direct 
method are beyond the capacity of present hardware. 

One way to circumvent the LBB condition is to modify the classical Galerkin mixed 
functional by appending some least-squares terms as Hughes and co-workers[12] did for 
their mixed Galerkin/least squares method. However, this type of method requires the 
choice of parameters, which depend on the mesh size of individual elements, and results in 
nonsymmetric matrices, which are still hard to deal with. An alternative approach to the 
use of equal-order elements (or unstaggered grids in a finite difference context) is to add 
the Laplacian of the pressure term and the divergence of the momentum into the continuity 
equation and to add a vector identity which relates the velocity and the vorticity into the 
momentum equations. Then standard finite volume, finite difference and Galerkin finite 
element methods can be employed to discretize the equations. This approach is proposed 
by Hafez and Soliman[13], and is essentially related to the mixed Galerkin/least-squares 
method. 

In recent years a great deal of attention has been drawn to projection methods 
which have numerous aliases: splitting methods, fractional step methods, pressure cor- 
rection methods, velocity correction methods and generalized and simplified marker-and- 
cell(GSMAC) methods. This family of methods was originated with the work of Chorin[14] 
in a finite difference setting, and was transferred to the finite element version, e.g., see 
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the studies by Gresho[15], Ramaswamy[16], Shimura and Kawahara[17], and Yagawa and 
Eguchi[18]. In these methods, the pressure computation is uncoupled with that of the 
velocity. The velocity is updated by explicit, implicit, or semi-implicit time-marching, and 
the pressure is obtained by solving the Poisson equation. These methods generally require 
less execution time and storage than the classical mixed Galerkin methods. This type of 
method should be classified as the mixed Galerkin method, and the LBB condition still 
should be respected. 


In the category of projection methods we should especially mention adaptive hp- 
methods developed by Oden and co-workers[19]. Their techniques can vary simultaneously 
the mesh size h and the spectral order p of elements to produce high resolution of flow 
features with a relatively small amount of computation. Their techniques combine some 
new general multistep schemes based on high-order characteristic methods with pressure 
projections which are handled through a combination of exterior penalty methods and an 
adaptive pressure correction technique. 


For the solution of incompressible Navier-Stokes equations, we have been developing 
a least-squares finite element method (LSFEM) [20-22], This method is based on the first- 
order velocity-pressure-vorticity formulation. Using C° finite elements to discretize the 
equations and minimizing the L 2 norm of the equation residuals lead to a symmetric and 
positive-definite algebraic system which can be effectively solved by simple matrix-free 
iterative methods. This is a minimization problem rather than a saddle point problem, 
and the choice of elements is thus not subject to the restriction of the LBB condition. In 
other words, all variables can be interpolated by the same element. This method is free 
of any parameters. There is neither added dissipation or upwinding, that is, the LSFEM 
is clean and robust. Since no derivatives are involved in boundary conditions, and only 
physical boundary conditions are imposed, the implementation of boundary conditions is 
extremely easy. In this method, all unknown variables (velocity, pressure and vorticity) 
are solved simultaneously, and there is no complicated iteration between them. 

The capabilities of the LSFEM have been further shown by Lefebvre et al [23] for 
two-dimensional unstructured triangular meshes, and by Tang and Tsang [24] for two- 
dimensional time-dependent flows with thermal convection. 


In the present paper the LSFEM is extended to large scale/ three- dimensional com- 
putation. Following this introduction, the velocity-pressure-vorticity formulation for in- 
compressible Navier-Stokes problems is presented in Section 2. An alternative velocity- 
Bernoulli function-vorticity formulation is discussed in Section 3. The numerical proce- 
dures for solving the first-order equations are addressed in Section 4. The numerical results 
for 2D and 3D cavity flows are contained in Section 5. Concluding remarks are given in 
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Section 6. 


2. Velocity-Pressure- Vorticity Formulation 

Let us consider the following steady-state incompressible Navier-Stokes problem in a 
bounded 2D or 3D domain f 1: Find the velocity u = (u,v,w) and the pressure p such that 


<1 

SI 

II 

o 

in fi, 

(la) 

u • Vu + Vp - ^ An = f 

in f h 

(16) 


Here all variables are nondimensionalized, f — {fnfyifz) is the body force, and Re denotes 
the Reynolds number, defined as 


v 


where L is a reference length, U a reference velocity and v the kinematic viscosity. Of 
course, boundary conditions should be specified to complete the definition of the boundary 
value problem. 

Since the momentum equations(lb) involve the second-order derivatives of velocity, 
direct application of the least-squares method requires the use of inconvenient C 1 elements 
and produces matrices with large condition number. In order to use the LSFEM with C° 
elements, we have to consider the governing equations of incompressible flow in the form 
of a first-order system. Therefore, we introduce the vorticity u> = (u> x ,w y ,u> z ) = V x u as 
an independent unknown vector, and rewrite the incompressible Navier-Stokes equations 
in the following first-order quasi-linear velocity-pressure-vorticity formulation: 


V • u = 0, 


(2a) 


u • V-u + Vp + — V x w = /, 


oj — V xu = 0. 


(26) 

(2c) 


The first-order system (2a), (2b) and (2c) has an odd number (seven) of unknowns, 
i.e. ( u,v,w,p , Wj.jWy.Wj), and an odd number (seven) of equations, i.e. one continuity 
equation, three momentum equations, and three definitions of the vorticity components. 
When we classify this first-order system( assuming that the convective terms have already 
been linearized) , we calculate the eigenvalues of a corresponding 7x7 coefficient matrix. 
Obviously, this matrix has at least one real eigenvalue which means that this first-order 
system cannot be elliptic in the ordinary sense. Therefore, the least-squares finite element 
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method based on (2a), (2b) and (2c) cannot have an optimal rate of convergence for all 
unknown variables. In order to avoid this trouble, we add the following compatibility 
condition into the system: 

V - w = 0, (2d) 

so that we have seven unknowns and eight equations, and the ellipticity can be proved by 
using the same technique as discussed in [25]. 

For Cartesian co-ordinates, the system (2) is given as 

du dv dw _ 

dx dy dz ’ 


du du du dp 1 . du z 

dx dy dz dx °- v 


_ ®fjL) ~ f 
Re y dy dz } Je ’ 


dv dv dv dp 1 , du x du z . 

“ fc + * Ty + W Tz + Ty + Te ( -§T - &■> = 


dw dw dw dp 1 .du y 

“ -di +V -d^ +W -dI + Tz + Te { ^- 

dv dw 

U * + Tz~ a? -0, 

dw du 

w ’ + d^-Tz =(> ' 


du x 

dy 


)»/« 


w z + 

du x 


du dv 
dy dx ’ 


dx dy 


du y duz _ Q 


dz 


(3) 


Since the system is of first-order, the boundary conditions are thus very simple, and do 
not involve the derivatives of unknowns. Let (ri,r 2 ,r 3 ,r 4 ,rs) denote the pieces of the 
boundary surface T. The unit outward normal vector to T is denoted by n, and the 
tangential vectors to T by Tj and r 2 - We may consider, for instance, the following boundary 
conditions: 

(a) u = 0,u = 0, w — 0 on Ti (the wall); 

(b) u = constant, v = 0,u> = 0 onr 2 (the far field); 

(c) u = given, v = 0 ,w = 0 on T 3 (the well developed inflow or outflow); 

(d) u T i = 0, u T 2 = 0 ,p = constant on Ts (the outflow); 

(e) u n = 0 ,u T i = 0, u T 2 = 0 ,p = constant on I\j (the free surface). 
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We note that in most cases the specification of boundary conditions for the vorticity 
components is not necessary. At sohd wall and well developed inflow or outflow boundaries, 
we prescribe only the velocity components, no vorticity is involved. 


3. Velocity-Bernoulli Function- Vorticity Formulation 

The system (2) is not the only way to write a first-order formulation by using the 
vorticity as an independent unknown vector [26]. By introducing the Bernoulli function(the 
total enegy) 6 = p + \(u\ 4- u% + u%) as an independent variable instead of the pressure p, 
we have the following Navier-Stokes problem: Find u, b and u> in Cl such that 


V • u = 0, 

(4a) 

-u x u> + V6 - V xuj — f, 

Re 

(46) 

ZJ — V x u — 0, 

(4c) 

o 

II 

13 

t> 

(Ad) 

This first-order velocity-Bernoulli function-vorticity formulation is 

also suitable for the 


LSFEM. 

For Cartesian co-ordinates, the system (4) is given as 

du dv dw_ 

Vx + d~y + ^~°^ 


db 


1 , du> z du)y 


+ di + Tc { ^~~dr ] - fx ' 

db 1 ,dw x dw z 

UUz - w. + Ty + fv. 


db 
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-UU, + - 

+ Ye 
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dv 

dw 

U? x + 

dz 

dy 

Vy + 

dw 

du 

dx 

dz 


du 

dv 

U> z + 

dy 

dx 


du> 2 
dy 


■) = /*. 


( 5 ) 
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dto x du) y dw z 

dx + dy dz 


For two-dimensional problems, let u? = w z , we have 


du dv _ 

9x 9y ’ 


96 1 9a ; 

+ a- + ^“'^T = /*> 

9x ite ay 
96 1 9w 

UU; dy Re dx y ' 


du dv 

" + a; ~ = °- 


( 6 ) 


The velocity-Bernoulli function- vorticity formulation has some advantages from a the- 
oretical point of view. In this formulation, the nonlinear terms u x u> are of zero-order, 
that is, they are not related to any derivatives, and the rest of the terms constitute a 
linear Stokes problem. Therefore, the whole system is elliptic. This fact may simplify the 
nonlinear analysis of the Navier-Stokes equations. 


4. Algorithm of Solution 

The quasi- linear problem (3) (or (5), or (6)) can be linearized by the successive substi- 
tution scheme, or the Newton’s scheme. For example, for the equations (3), the nonlinear 
convection term can be linearized as u° by using the simple substitution, or as 

u° -f — u °^- by using the Newton linearization scheme. Here the superscript “0” 
indicates that the value of the corresponding variable is taken from the previous calculation 
step. For the equations (6), the nonlinear term — vw can be linearized as — v°w by using 
the simple substitution, or as -v°w - vw° + v°iu 0 by using the Newton’s scheme. 

The linearized first-order equations are now treated by the LSFEM[20]. The LSFEM 
results in symmetric and positive-definite algebraic equations. A Jacobi preconditioned 
conjugate gradient method is then employed to solve the linear algebraic equations. In 
the conjugate gradient method, the major computation is the multiplication of the global 
matrix with the global vector, and this can be done in an element-by-element manner 
without forming the global matrix[27]. In order to further save the storage, in our algorithm 
even the element matrices are not formed, We directly calculate the product of the element 
matrix and the element vector. At the same time, the Jacobi preconditioner, which consists 
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of diagonal terms of the global matrix, can be easily formed. In this way we store only 
several global vectors, and the derivatives of the shape functions at Gaussian points for 
each element. 


5. Numerical Examples 

5.1. Comparison of Two Formulations 

To compare the performance of the velocity-pressure-vorticity formulation and the 
velocity- Bernoulli function-vorticity formulation, we performed a series of computations 
for 2D cavity flow at different Reynolds numbers with 50 x 50 non-uniform bilinear ele- 
ments(the mesh is the same as that used in [21]) using the following four schemes: 

(a) Scheme PS, the (u,v,p,o;) formulation with simple substitution; 

(b) Scheme PN, the (u, v, p, w) formulation with Newton linearization; 

(c) Scheme BS, the (u,v,b,uj) formulation with simple substitution; 

(d) Scheme BN, the (u,v,6,u>) formulation with Newton linearization. 

In order to eliminate any effects from iterative solvers, in all these comparison studies 
a direct solver was employed to solve the resulting linear algebraic system. Our numerical 
results reveal that, as for the speed of convergence, Scheme PS and Scheme BS are almost 
comparable, and Scheme PN and Scheme BN are almost comparable. However, in some 
cases, the schemes based on the (u,v,b,u) formulation need a little more iterations or a 
closer initial guess to converge. Figure 1 shows the convergence history of the I^-norm of 
residuals for the cavity flow at Re = 5000 using different schemes. For Scheme PS and PN, 
the initial guess of u and v is taken from the results of Re = 3200. For Scheme BN the 
computation started from the results of Re — 4000. If the results of Re — 3200 are used as 
the initial guess, Scheme BN does not converge. Here Scheme PS was implemented with a 
relaxation number of 0.8[21] . IVithout relaxation Scheme PS does not converge. Scheme 
BS does not converge even with a relaxation number of 0.8 and the use of the results at 
Re = 4000 as the initial guess. Based on the results of our numerical experiments, we are in 
favor of the (u,v,p,u>) formulation. In the following computation of large-scale examples, 
Scheme PS is used with the Jacobi preconditioned conjugate gradient method. 

5.2. 2D Driven Cavity Flow 

The LSFEM solution of the 2D cavity problem has already been reported in [21] by 
using the steady-state approach and in [24] by a time-dependent algorithm. In these pre- 
vious studies, a grid of 50 x 50 nonuniform bilinear elements was employed. The numerical 
results at Re = 10000 are in good agreement with the fine mesh 257 x 257 results of Ghia 
et al.[28]. 

In this work, the objectives of choosing the 2D driven cavity problem are twofold. 
First, we would like to show the capability of the LSFEM for large-scale problems. Second, 


8 


we would like to find out whether improved results can be obtained by using very fine grids. 

The definition of driven cavity flow is as usual. The boundary conditions for (u,v) 
are u = v = 0 on all solid walls and u = l,v = 0 on the top lid. Also p = 0 is specified 
at the centre of the bottom. No boundary conditions are given for the vorticity. At first 
we divided the cavity into 50 x 50 uniform bilinear elements with the size h = 0.02. In 
order to take care of the singularity at the top corners, we again divided the top layer of 
elements into two thin layers of thickness 0.005 and 0.015 units, respectively. 

At first we calculated the flow at Re = 10000 with this 50 x 51 mesh. By using a kind 
of continuation method discussed in [21], the computation could be done on a IBM PC-386 
with 4M bytes memory in one day. Then we refined the mesh to 100 x 102 elements, and 
interpolated the previous solution as the initial guess for the conjugate gradient method. 
We applied this procedure sequentially until the steady-state solution was obtained for the 
400 x 408 mesh in which most element have the size h = 0.0025, and in the top layer the 
thickness of the elements is 0.0005. The storage required for this problem is less than 8M 
words on CRAY-YMP. The residual of each discretized equation at Gaussian points was 
less than 10 -6 . 

The profile of the horizontal velocity for Re = 10000 along the vertical center line 
of the driven cavity (x = 0.5) is illustrated in Figure 2(a). Overall, the profile compares 
well with that given in Ghia et al. [28], except in the boundary layers. In our results, the 
boundary layer phenomena are more pronounced. The computed streamlines and vorticity 
contours are given in Figure (2b) and (2c). 

5.3. 3D Driven Cavity Flow 

We choose the cubic driven cavity problem shown in Figure 3 to further test our 
method. This problem has been simulated by using finite difference[13, 29-36] and spec- 
tral methods[37]. Due to the difficulties of existing finite element methods for large-scale 
computation, it may not be surprising that finite element solutions for this problem are 
scarce. Only Kato, Kawai and Tanahashi [38] reported the detail of numerical results for 
cubic cavity flows by using their GSMAC finite element methods. It should be noted that 
Kato et al. use the Bernoulli function and the vorticity to write the momentum equations. 
However, in their computation the vorticity vector is not independent, but calculated from 
the velocity field. 

Most of numerical efforts were restrictive in mesh size and insufficient in resolution, 
due primarily to the limitations on numerical methods and computer resources, except 
Iwatsu et a/.[35,36], who employed very fine stretched mesh(81 x 81 x 81). 

Almost all of the previous researchers used time-marching schemes, and claimed that 
steady-state solutions were obtained by their own standards of convergence for Re < 2000, 
except Rosenfeld et a/. [32] mentioned that their solution at Re = 100, 1000 was still varying 
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in the lower corners. Also no Taylor-Gortler-like vortices were reported for Re < 2000 by 
previous researchers. 

It should be mentioned that Koseff and Street [39-41] conducted systematic experi- 
ments of flow in a driven cavity. Their experiments revealed that the flow field is highly 
unsteady and possesses significant secondary motions (end-wall corner eddies and Taylor- 
Gortler-like vortices). However, in all of their papers, they showed only the results for 
Re = 3200, and a very important observation, that is, for what Reynolds number the flow 
becomes unstable, had never been reported. 

In the present study, we carried out computation at Re = 100,400,1000,2000, and 
3200. u = 1, v = to = 0 were specified on the top driven surface ( y = 1.0), and u = v = 
w = 0 on all solid walls(Figure 3). At the center of the bottom (* = 0.5, y = 0.0, z = 0.5), 
p = 0 was specified. There were no boundary conditions on the vorticity. Three meshes 
were used. The first mesh consisted of 50 x 50 x 50 non-uniform trilinear elements. The 
smallest element size had h = 0.002 at the corners; the largest had h = 0.04 at the center. 
The second mesh had 50 x 50 x 50 uniform trilinear elements. The third mesh was based 
on the second one. In order to take care of corner singularities two layers of thin elements 
were added close to the top driven surface. So the third mesh had 50 x 52 x 50 elements. 
The problems were solved using less than 13M words of memory on a CRAY-YMP. 

The observations from our numerical results can be summarized as follows: 

(a) Up to Re = 3200 the flow is symmetrical about the plane z = 0.5. This observation 
is consistent with the results of numerical simulation published by other researchers. 

(b) The 3D cavity flow is highly complicated. For 2D driven cavity problems the small 
eddies appear along the boundaries. It is a common practice to use fine grids along the 
boundaries to capture the flow details. However, for 3D problems a pair of small vortices 
may be formed near the center part of the cube even for flows with moderate Reynolds 
numbers (Re > 1000). Therefore, stretched 3D meshes may not be suitable for high- Re 
flows. It means that for 3D simulation we should use uniform fine grids to capture small 
eddies both inside and outside of the boundary layer. This makes 3D computation much 
more difficult than 2D computation. 

(c) The 3D cavity flow is highly unstable. For 2D flows at high Re number(10000 
for the cavity flow) by using the LSFEM it is quite easy to reduce the L 2 norm of the 
residual(assuming the whole area of the flow domain is 1) to the level of 10 -5 . However, 
for the 3D cavity flow, we were not able to obtain the steady-state solutions for Re > 1000. 

(d) The Taylor-Gortler-like vortices are observed for Re > 1000. 

Some numerical results are provided in Figure 4, 5 and 6. Figure 4 and 5 show the 
projections of velocity vectors and the contours of vorticity components on each section at 
Re = .100 and 400. The flow patterns at Re = 100 and 400 are in good agreement with 
the results of Iwatsu et al. [35,36]. The vorticity contours compare well with those of Hafez 
and 
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Soliman[13]. Figure 6 clearly shows the corner vortices and TGL vortices at the bottom 
region of the cavity at Re = 1000. 


6. Conclusions 

A least-squares finite element method based on the velocity-pressure-vorticity formu- 
lation is successfully extended to the solution of large-scale/3D steady-state incompress- 
ible Navier-Stokes problems. We prefer the velocity-pressure-vorticity formulation to the 
velocity-Bernoulli function- vorticity formulation, because the former is faster and not sen- 
sitive to initial guesses. The LSFEM generates a symmetric, positive-definite algebraic 
system of equations which can be robustly solved by the matrix-free conjugate gradient 
method. In this method there is neither upwinding, adjustable parameters, numerical 
boundary conditions, splitting, projection, nor artificial compressibility. Besides the finite 
element interpolation and the linearization, no other approximation is introduced into this 
method. Therefore, the solution of this method is more reliable than that of any existing 
method. 

By using this method the steady-state solution is computed for 2D cavity flow at 
Re = 10000 with 400 x 408 bilinear elements. Our results suggest that the benchmark 
solution of Ghia et al. [28] may not be accurate enough in the boundary layer. For 3D driven 
cavity flows the presence of Taylor-Gortler-like vortices is observed at Re > 1000. We hope 
to confirm these results with further works on the unsteady Navier-Stokes equations. 
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Figure 3- Geometry and boundary conditions of cubic cavity flow 





(d) Velocity vector at y=0.5 (e) VorticityK) contours at y=0.5 (f) Pressure contours at y=0.5 
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(a) Velocity vector at x=0.5 (b) Vorticity(o) x ) contours at x=0.5 (c) Pressure contours at x=0.5 



(d) Velocity vector at y=0.5 (e) Vorticity(w y ) contours at y=0.5 (f) Pressure contours at y=0.5 



(g) Velocity vector at z=0.5 (h) Vorticity(w 2 ) contours at z=0.5 (i) Pressure contours at z=0.5 


Figure 5. Numerical results for 3D cavity flow at Re = 400 
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(d) Velocity vector at y=0.5 (e) Vorticity(u> y ) contours at y=0.5 (f) Pressure contours at y— 0.5 



(g) Velocity vector at z=0.5 (h) Vorticity(a> z ) contours at z=0.5 (i) Pressure contours at z=0.5 

Figure 6. Numerical results for 3D cavity flow at Re = 1000 
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